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ABSTRACT 

In this paper we consider the evolution of small planetesimals (radii '^1 — 10 metres) 
\ in marginally stable, self-gravitating protoplanetary discs. The drag force between the 

. disc gas and the embedded planetesimals generally causes the planetesimals to drift 

inwards through the disc at a rate that depends on the particle size. In a marginally 
stable, self-gravitating disc, however, the planetesimals are significantly influenced by 
the non-axisymmetric spiral structures resulting from the growth of the gravitational 
' instability. The drag force now causes the planetesimals to drift towards the peaks 

, of the spiral arms where the density and pressure are highest. For small particles, 

that are strongly coupled to the disc gas, and for large particles, that have essentially 
decoupled from the disc gas, the effect is not particularly significant. Intermediate sized 
particles, which would generally have the largest radial drift rates, do, however, become 
O ' significantly concentrated at the peaks of the spiral arms. These high density regions 

^ \ may persist for, of order, an orbital period and may attain densities comparable to that 

C/3 i of the disc gas. Although at the end of the simulation only 25 % of the planetesimal 

' particles lie in regions of enhanced density, during the course of the simulation at 

least 75 % of the planetesimal particles have at some stage been in a such a region. 
, , We find that the concentration of particles in the spiral arms results in an increased 

' collision rate, an effect that could significantly accelerate planetesimal growth. The 

density enhancements may also be sufficient for the growth of planetesimals through 
direct gravitational collapse. The interaction between small planetesimals and self- 
gravitating spiral structures may therefore play an important role in the formation 
of large planetesimals that will ultimately coagulate to form terrestrial planets or the 
cores of gas/ice giant planets. 

Key vifords: accretion discs — stars:pre-main-sequence — planetary systems: pro- 
toplanetary discs — planetary systems: formation 



1 INTRODUCTION 

The formation of both terrestrial planets and the cores 
of gas/ice giant planets is thought to occur through the 
coUisional accumulation of planetesimal s (Safronov 1972; 
IWetherilj[l990l:IWeidenschiUing fc Cuzzi|[r993D . In the case 
of gas giant planets, a gaseous envelope is accreted once 
the core has become sufficiently massive llLissaueil Il993t 
IPoUack et al]|l99d) . Since this must occur while there is 
still e nough gas in the cir cumstellar disc, and since observa- 
tions jHaisch et al.ll200lh suggest that most stars lose their 
gaseous discs within ~ 10^ years, it is generally accepted 
that gas giant planets have to form within ~ 10^ years. 

The standard core accretion models of giant planet for- 



mation jPoUack et aljll996l : iBodenheimer et al.ll200ol) sug- 
gest formation times that could easily exceed disc lifetimes. 
However, this long formation timescale problem might be 
solved in light of the results of recent numerical simula- 
tions of the evolution o f planetary cores embedded in tur- 
bulent accretion discs jNelson fc Paoaloizoul 120041) . These 
simulations show that cores may actually undergo a ran- 
dom walk through the disc, leading to the suggestion that 
core migration may significantly accelerate gas g iant forma- 
tion jRice fc Armitagel2003l:lAUbert et alJl2004l). a possibil- 
ity initiallv recognised bv lHourig£^^^W^J^1984^ . 

An additional difficulty in the standard core accretion 
model is the growth of kilometre sized planetesimals from. 
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initially, micron sized dust grains. Within a few scaleheights 
of the disc midplane, the pressure gradient in the circumstel- 
lar disc tends, for standard disc geometries, to be negative 
and causes the gas to orbit with sub-Keplerian velocities. 
The dust, which is not affected by the gas pressure gradi- 
ent, and the gas therefore orbit with different velocities and 
the resulting drag force causes the dust gr ains to drift in- 
ward s at a rate th at depends on their size iWeidenschillinel 
Il977l: l Takeuchi fc Lin 2002il . For small sizes, the dust grains 
are essentially coupled to the disc gas and the radial drift 
velocity is consequently small. For large sizes, the grains 
are decoupled from the gas, move in nearly Keplerian or- 
bits, and again have small radial drift velocities. Particles 
with intermediate sizes can, however, have large inward ra- 
dial velocities. Although the exact size range depends on the 
circumstellar disc properties, the maximum radial velocity 
may easily exceed lO'^cm/s and is normally thought to occur 
for ob jects with sizes between 1cm and Im llWeidenschillinel 
119771) . 

Although the differential radial velocity is a cru- 
cial part of the grain growth, since larg er o bjects grow 
by sweeping up smaller ob jects (e.g., ISafro nov 1 973 
IWeidenschiUing fc Cuzzilll993D . if the maximum radial ve- 
locity is too high these objects may drift inwards before 
they can become large enough to decouple from the disc 
gas. Together with the inward radial dri ft, the l arger grains 
may also settle towards t he midplane llGoldreich fc Ward 
Il973l: iGaraud et aT]l2004) . producing a thin, dense dust 
layer (note, however, that the presence of t urbulence in 
the disc can prevent the dust settling, see IStone et alJ 
Il99d) . This layer may become self-gravitating and, if suf- 
ficiently unstable, could produce kilometre sized objects d i- 
rectly via gravitational collapse jGoldreich fc Wardl[l973) . 
This, however, requires extremely small random dust ve- 
locities (^ 10 cm/s), which may be difficult t o achieve 
JWeidenschiUing fc Cuzzilll993l : ICuzzi et alJll993l) . Even if 
such small random velocities are possible, if objects are to 
grow via gravitational collapse, any incre ase in the ran- 
dom velocities has to be lost very rapidly llGammid 120011: 
[Rice ct al. 2003a). 

An alternative mechanism for the formation of gas giant 
planets, that does not require the initial growth of a core, is 
that the gas itself may become unsta ble, producing gir avita- 
tionally bound gaseous protoplanets iBossll998ll200(](l . This 
again requires th at the gas be cold and th at it remain 'al- 
most isothermal ' JPickett et alJll998l . l2000D or, equ ivalently 
that any cooling mechanis m is extrem ely efficient iGammij 
[2OOI; Rice et al. 2003a; Joh nson fc G ammic 20q3). Further- 
more, some simulations of fragmenting protoplanetary discs 
suggest that, at best, this mechanism may be able to pro- 
duce onl y the most massiv e (> 5 Jupiter masses) gas giant 
planets faice et al.ll2003b^ . 

We study here the influence of the development of grav- 
itationally unstable spiral modes on the planet formation 
mechanism via core accretion (see also lHaehighipour fc Bosa 
liooli). Although it is possible that the condi tions required 
for d isc fragmentation may never be met jPickett et alJ 
I2OOOI) . it is quite likely that protost ellar discs are self- 
gravitating in their ea rly stages (e.g.. iLin fc Pringl3ll990t 
iLodato fc Bertinll200lD . If so, the disc then becomes suscep- 
tible to the growth of non-axisymmetric spiral structures 
which can transport angular momentum very efficiently 



(e.g., iLodato fc rUc3 l2004l) . In the presence of such spiral 
structures, the gas pressure gradient, which can be large, 
changes from positive on one side of the structure to negative 
on the other, resulting in both super- and sub-Keplerian gas 
velocities. The gas drag force then causes dust grains to drift 
both radially inwards and outwards, depending on whether 
the local gas velocity is supe r-Keplerian or sub-Keplerian 
jHaghighipour fc Bosj|2003allbl) . The net effect is that the 
dust drifts towa rds the density maxima, wher e the pressure 
gradient is zero llHaghighipour fc Bosj^2003a^ . A similar ef- 
fect would occur in the presence of any coherent an d long- 
lived density enhancement. Previo us work (Godon fc Livid 
I 2000l : lKlahr fc Bodenheimeill2003l) has considered how vor- 
ticies may influence embedded particles. The formation of 
such vorticies is, however, still a matter of some debate. 

In this paper we present the results of three- 
dimensional, global simulations of self-gravitating accretion 
discs in which we include bo t h gas and dust, coupl ed via 
a drag force JWhippld IT97i: IWeidenschiUind Il977ll . The 
gaseous disc is maintained in a state of marginal gravita- 
tional instability, achieved by letting the disc cool down 
(through a sim ple pajametrisation o f the cooling function; 
ICammielboOll and lRice et "31112003311 . In this way a quasi- 
steady spiral structure develops in the disc and is maintained 
during several dynamical time-scales. We are therefore able 
to follow the process of concentration of the planetesimals 
in the spiral arms. We have performed several simulations, 
considering planetesimals of different sizes. We find that, 
for a given size range, the planetesimals are indeed able to 
reach high concentrations, in some regions attaining densi- 
ties comparable to the gas. This could significantly enhance 
the coagulation of planetesimals into larger bodies, by in- 
creasing the planetesimal collision rates and/or by making 
the planetesimal sub -disc be come gravitationally unstable 
( Youdin & Shu 200i lYoudirr fc Chiang 20o4- 

This paper is organised as follows. In Section|21we sum- 
marize the basic features of the coupled dynamics of a two 
component (gas and planetesimals) disc, including a descrip- 
tion of the relevant drag forces. In Section|3]we describe the 
simulation code that we have used and the numerical setup. 
In Section^ we describe our results. In Section]^ we draw 
our conclusions. 



2 GAS-PLANETESIMAL DISC DYNAMICS 

We consider a system comprising two interpenetrating discs: 
a "gas" disc, that is evolved using the standard hydrodynam- 
ical equations of motion, and a "planetesimal" disc, that is 
considered as a collection of test particles evolved under the 
influence of gravitational and drag forces alone. Both the 
"gas" disc and the "planetesimal" disc rotate around a cen- 
tral protostar of mass A/*, which we take to be equal to 
IMq. The term "planetesimal" is generally used to refer to 
particles that have decoupled from the disc gas. Although 
the particles we consider here are still coupled to the gas, 
we use this term because, as we will discuss later, the parti- 
cles we consider are clearly too large to be regarded as dust 
grains. 

In order to illustrate the basic dynamical ingredients of 
our model, and to introduce the relevant physical quantities, 
let us consider first the simple case of a smooth, axisym- 
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metric, non self-gravitating disc. We define the Keplerian 
velocity, Vk, around the star at a distance R by 

Vl = GAL/R. (1) 

The gas disc is characterized by a surface density profile 
S(_R), a temperature profile T{R) and a pressure profile 
P{R). The gas sound speed is defined as 



dP 
dp' 



(2) 



where p is the gas volume density. Let v and Vp be the 
velocities of the gas and of the planetesimals, respectively. 
In centrifugal equilibrium the azimuthal component of the 
gas velocity is given by 



2 

V4, 



d In p 
dlni? 



^i(i + 7), 



(3) 



where 7 = (dlnp/dln J?)Cs/V^, is a measure of the impor- 
tance of thermal effects in the disc. For most disc models, 
the density p decreases with radius, so that 7 < and «^ is 
generally sub-Keplerian. If the gas disc is in vertical hydro- 
static equilibrium, its aspect ratio is given by 
-1/2 

VW\ (4) 



'r 



Cs 

Vk 



d In p 



dlni? 



The planetesimals, on the other hand, are not affected 
by pressure forces and, in centrifugal equilibrium therefore 
orbit at the Keplerian velocity Vk- Let u be the relative ve- 
locity between the gas and the planetesimals. We assume 
that the gas and the planete simals are coup l ed thr ough 
a drag force, a s described by Wcid enschillind ^1977^ and 
IWhiPDlel (|l972l ). The drag force is given by 

Fd = —^Cvna^pu^u, (5) 

where u = |u|, u = u/m, a is the mean radius of the plan- 
etesimals and Cd is the drag coefficient, given by 



Cd = < 



3 u 
24R-^ 

0.44 



a < 9A/4 

Re < 1 
l<Re< 800 
Re > 800 



(6) 



The drag regime where a < 9 A/4 is generally called the 
Epstein regime, whereas the other three regimes define the 
Stokes drag. In the previous expression. Re is the Reynolds 
number, defined below in Eq. ||5J and A is the mean free path 
of the gas particles. Assuming the gas to be made mainly of 
molecular hydrogen, this is given by 



A = 



4 10" 



(7) 



pA p 

where in the last expression the gas density p is evaluated 
in cgs units, uih-, is the mass of the hydrogen molecule, and 
A is its cross section 



A = TTal^7 10-^'^cm^ 



(8) 



where gp is the mean radius of the hydrogen molecule (see 
ISupulver fcTinlbOOOTl . 

The Reynolds number Re is given by 



Re — 



2apu 



(9) 



where rj = pi^ is the gas viscosity. For coUisional viscosity, 
we have 
pcsA 



We can therefore write 



Re 



(10) 



(11) 



For the disc properties and particle sizes that we will con- 
sider here, the Mach number rarely exceeds unity and the 
Reynolds number generally falls in the range 1 < _Re < 800. 
Only in the inner regions of the disc, where u is largest, 
does the Mach number exceed unity and the Reynolds num- 
ber exceed 800. 

Another important quantity is the "stopping time" te, 
where 



tc = 



Fd 



(12) 



Here nip = 47rpsa^/3, is the mass of the planetesimals and ps 
is the internal density of the planetesimals (Throughout the 
paper, we take ps = 3 g/cm"^). With the previous definitions, 
we have 



8 (Ps 
3Cd V p 



(13) 



where td ~ f^K^ ~ R/Vk is the local dynamical timescale. 
The smaller ta/t^, the more strongly the planetesimals are 
coupled to the gas. It is well known j Weidenschillins| 1 1 9 77^ 
that the gas-planetesimal interaction through the drag force 
causes the planetesimals to migrate in the direction of in- 
creasing pressure (i.e., inward, for a smooth, axisymmetric 
disc). The migration rate depends on the particle size and 
there is a range of particles sizes (the extent of which de- 
pends on the gas disc properties) that are characterized by 
a relatively large radial drift. 

Consider, as an illustration, the average structure of one 
of the gas discs obtaine d from the num e rical simulations of 
self-gravitating discs bv lLodato fc Ric3 ((ioO^). On average, 
these discs are characterized by a power-law surface density 
profile, E oc R~^, and by an approximately flat profile, with 
a value of order unit y, of the axisy mmetric gravitational sta- 



a value ot order unit y, 01 tne axisy ir 
bUity parameter Q iToomrelll964l') 



(14) 



In the previous expression k is the epicyclic frequency, that, 
for the nearly Keplerian discs considered here, is roughly 
equal to the orbital frequency Q,, to order H^/R^. The disc 
extends from Rin ~ 0.25 an to i?out ~ 25 au. The total 
disc mass is Mdisc = O.25M0. The profile of Q essentially 
determines the pressure structure in the disc. We have then 
that Cs cx R^^'^, H/R oc R, and p cx R~^. 

We can therefore estimate the radial drift velocity of 

S ilanetesimals in such a disc, using standard techniques 
Weidenschilline^ll977^ . The results are shown in Fig. 
where the left panel shows the radial drift velocity Ur in 
units of the local Keplerian velocity. Equivalently, defining 
the drift timescale U = R/ur, we have Mr/Vk = t^/tr, so 
that the left panel of Fig. also shows a measure of the 
drift timescale compared to the dynamical timescale. In the 
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Figure 1. Left panel: Radial drift velocity (in units of the local Keplerian velocity) as a function of planetesimal size for the assumed 
disc structure. Right panel: Radial drift timescale as a function of planetesimal size. The curves refer to four different locations in the 
disc: (solid line) 1 au, (dotted line) 3 au, (dashed line) 5 au, (long-dashed line) 10 au. 



right panel of Fig. we also show ty in yrs. It can be seen 
that for the assumed average gas disc properties, the largest 
radial drift would occur for planetesimals with sizes in the 
range between 10 and 100 cm. It should however be kept 
in mind that these results refer to the azimuthally averaged 
structure of the disc. Actually, the self-gravitating gas disc is 
characterized by a spiral structure, with alternating regions 
of high and low gas densities (see Fig.|21below). Equation ^ 
then readily shows that the gas velocity changes from sub- 
Keplerian through super-Keplerian, when moving across one 
arm of the spiral structure. The drag force causes the plan- 
etesimals to drift towards regions of higher pressure (and 
hence of higher density), so we expect that the planetesi- 
mals tend to concentrat e very fast around the maxi ma of 
the gas density (see also iHaghighipour fc Bosjl2003a|) . The 
results described above then suggest, for the disc we are con- 
sidering here, that this effect is maximal for planetesimals 
with sizes between 10 and 100 cm. 

The expected drift timescale, in such structured discs, 
is going to be smaller than in a smooth disc. In fact, in a 
smooth, power-law disc the pressure gradient VP ~ P/R, 
while in a disc with a spiral structure (characterized by 
a typical scale ~ H) the pressure gradient is of order 
VP ~ P/H. The increased pressure gradient decreases the 
drift timescale by a factor H/R. The drift timescale is de- 
creased by further factor H/R because the planetesimals 
have to move only a distance of order H to reach the den- 
sity maxima. We therefore expect a significant concentration 
of the planetesimals around the gas density maxima to take 
place over a timescale a factor (H/R)^ ~ 0.01 smaller than 
that displayed in Fig.Q becoming therefore comparable with 
the dynamical timescale in the most favourable cases. By in- 
spection of the right panel of Fig. we also obtain that for 
particles with sizes between ~ 10 and ~ 1000 cm, we need 
only run our simulation for about 100 yrs in order to follow 
the effect of radial drift on the planetesimal sub-disc. 



3 NUMERICAL SIMULATIONS 

3.1 Smoothed particle hydrodynamics 

The three-dimensional gaseous disc used in these simula- 
tions is modelled using smoothed particle hyd rodynamics 
(SPH), a L agrangian hydrodynamics code (see lBenj[l99ol : 
iMonaghant 19921. The gas disc consists of 250,000 SPH parti- 
cles, each of which has a mass and an internal energy (tem- 
perature). Each particle also has a smoothing length that 
is allowed to vary with time to ensure that the number of 
neighbours (SPH particles within 2 smoothing lengths) re- 
mains ~ 50. These neighbouring particles are used to deter- 
mine the density which, together with the internal energy, is 
used to compute the pressure. The central star is modelled 
as a point mass onto which gas particles may accrete if the y 
approach to within the sink radius (e.g.. lBate et aljll995^ ■ 
here taken to be 0.25 au. Both the central point mass and 
the gas particles use a tree to determine gravitational forces, 
and to determine the gas particle neighbours. 

The small planetesimals are modelled using an addi- 
tional type of particle. These particles are, as far as the gas 
simulation is concerned, massless (i.e., in these simulations 
we neglect the self-gravity of the planetesimal disc and the 
back reaction of the drag force on the gas). They experience 
only gravitational forces (from the central star and from the 
disc gas) and are coupled to the disc gas via a drag force. 
As discussed in Section |5| the drag force depends on the 
particle size, the local gas density, and on the local gas ve- 
locity. To determine the drag force coefficient also requires 
the gas sound speed which is calculated using the local gas 
internal energy. The gravitational force on these test par- 
ticles is computed by including them in the tree. This also 
determines their nearest gas neighbours. To ensure that the 
number of gas neighbours remains ~ 50, the test particles 
also have a smoothing length that is allowed to vary with 
time. These neighbouring particles are then used to calcu- 
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late the gas density, velocity, and internal energy at the loca- 
tion of_each_test particle using the standard SPH formalism 
(see iM^niihiBEiii) • 

The exact value of the drag force is 
then determined by specifying the planetesimal size. In each 
simulation performed here, we use 125,000 test particles to 
represent a planetesimal disc which we assume contains par- 
ticles of a single size. 

An additional saving in computation al time is ma de by 
using individual particle time-steps (Ba te et alJIlQAsI ) with 
the time-steps for each particle li mited by the Cou rant con- 
dition and by a force condition llMonaghanlll993) . As will 
be discussed in more detail in the following sections, the gas 
disc is assumed, in the absence of any cooling mechanism, to 
have an adiabatic equation of state with an adiabatic index 
of 7 = 5/3. The gas disc is initially evolved, in the absence 
of test particles, by imposing a cooling term which is chosen 
such that the disc ultimately settles into a quasi-steady, self- 
gravitating state. The test particles are then added and the 
simulation is evolved for an additional outer rotation period. 



3.2 Self-gravitating gas disc simulations 



In this w ork we use the results of the numerical simu- 
lations by iLodato fc Rlcel ()2004 ) of the dynamics of self- 
gravitating gas discs as an input for our two-component 
gas-planetesimals disc simulations. 

It is well known that the onset of gravitational instabil- 
ities in the disc is determined by the value of the parameter 
Q, defined in Eq. I|14|l . If Q is smaller than a threshold value 
of order unity, the disc quickly develops a spiral structure 
on the dynamical timescale. The presence of the spiral in- 
fluences strongly the thermal evolution of the disc, in that 
it provides a source of effective heating. Since Q is propor- 
tional to the thermal speed Cs, in the absence of any cool- 
ing Q would rapidly become relatively large and the spiral 
structure would vanish. However, if some cooling is present, 
a self-regulated state can be achieved where the heating pro- 
vided by the spiral structure balances the external cooling, 
leading to a long-lasting spiral. 

Th ese processes have be en recent ly explored nu meri- 
cally bv lLodato fc Rlcel lj2004h (see alsolSimmiJiQQ^), who 
have performed global, three-dimensional simulations of the 
evolution of self-gravitating discs. They imposed a cooling 
of the form: 



dt/ 
dt 



U 



(15) 



where U is the internal energy of the gas and the cooling 
timescale is taken to be simply proportional to the dy- 
namical timescale, tcooi ~ For very rapid cooling 
timescales (/3 S, 3) the sel f -gravitating disc u ndergoes frag- 
mentation (|Gammiell200ll : iR.ice et alJl2003a^ . In this work 
we have taken /? = 7.5, giving a cooling time that should 
not, and indeed does no t, lead to fragmention . 

The simulations bv lLodato fc Ric3 (|22q3) indeed show 
the effectiveness of the self-regulation process. The disc ex- 
tends from 0.25 au to 25 au, and is characterized initially by 
a surface density profile E oc and a temperature profile 
T oc R^^^^. The exact surface density is determined by spec- 
ifying a total disc mass, and the temperature is determined 
by specifying that the Toomre Q parameter has an initial 
value of 2 at the outer edge of the disc. The temperature 



Figure 2. Surface density structure of a self-gravitating disc with 
A^disc = O.25M0 (the central star mass is M* = IMq). The 
spiral structure is a quasi-steady feature lasting for at least several 
thermal timescales. The panel shows the logarithm of the surface 
density S with the scale covering 1 < log(S/g cm"'^) < 4.7. The 
size of the box is 50 au across. 



profile, however, is rapidly modified by the competing heat- 
ing and cooling processes operating in the disc, as discussed 
above. At the end of the simulations a self-regulated state is 
achieved with an almost constant profile of Q, with a value 
close to unity. The spiral structure obtained in this way is 
a quasi-steady feature lasting for at least several thermal 
timescales (i.e., at least until the end of the simulations). 
Fig. |5| shows the final disc structure for a disc whose total 
mass is Mdisc = O.25M0. The image shows the logarithm of 
E, with a colour scale covering 1 < log(E/g cm~ ) < 4.7. 

The spiral structure transports angular momentum in 
the disc and therefore promotes the accretion process. This 
leads also, in some cases, to a steepening of the surface 
density profile. This process however occurs on the much 
longer 'viscous' timescale, so that the final profile of E is 
only slightly modifie d with respect to the initial one (see 
iLodato fc Ric^l2004^ . 



3.3 Numerical setup 

We initially evolve the gas disc in the absence of any test 
particles, as described in the previous Section. After ~ 6 
orbital periods at the outer edge of the disc (i.e., after ~ 800 
yrs), a quasi-steady, self-regulated state is reached. We then 
introduce the planetesimals and evolve the simulation as 
described above in Section 13.11 for roughly one more outer 
orbital period (i.e., ~ 125 yrs). 

In order to cover a wide range of planetesimal sizes we 
have considered separately planetesimal sizes of 50 cm and 
1000 cm. Based on the results of Section |21 we expect the 
50 cm sized planetesimal to have the largest radial drift. For 
smaller planetesimal sizes (e.g, 1 cm) we expect the accel- 
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Figure 3. Surface density structure of the planetesimal discs after one outer orbital time. The upper panels show Sp, for the 50 cm 
planetesimals (left panel) and for the 1000 cm planetesimals (right panel). The surface densities have been multiplied by 100 and the 
colour scaling is the same as in Fig. |^ in order to have a direct comparison with the gas surface density structure. The size of the 
box displayed is also the same as in Fig. \^ The bottom panels show the ratio Sp/S for the two cases (the colour scales of these 
panel are the same in the two cases and cover 0.004 < Sp/E < 0.04). If the gas and planetesimals respond in the same way to the 
gravitational instabilities, the density ratio would be uniform through the disc, so that any non-uniformity in the bottom panels is a 
direct measure of the concentration effect on the planetesimals caused by the combination of gas gravity and drag force. The planetesimal 
tend to concentrate in clumpy regions along the spiral arms, where the minima of the gravitational potential are located. The effect is 
particularly evident for planetesimal size of 50 cm, for which Sp/S can be enhanced by more than a factor of 50. 



eration due to the drag force to be very large, resulting in 
extremely long computation times. However, for such small 
planetesimals, we expect the drag to be so large that their 
structure will closely match that of the gas. In addition to 
the 50 and the 1000 cm cases, we have also performed one 
simulation in which no drag force was included, so as to 



provide a direct measure of the effect of gas drag on the 
evolution of the planetesimals. 

The planetesimal disc initially extends from i? = 2 au to 
R = 2Q au. The surface density profile of the planetesimals 
Ep was taken to be proportional to . Since we neglect 
the planetesimals' self-gravity and the back reaction of the 
drag force on the gas, the actual value of the planetesimal 
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Figure 4. Left panel: volume density distribution of the planetesimal for the two different sizes considered (solid line; 50 cm; dotted 
line: 1000 cm), at the end of the simulation. Note that the particular shape of the density distribution is not very meaningful, since it 
depends on the initial surface density distribution that we assume. What is interesting is that, starting from the same initial density 
distribution, at the end of the simulation the 50 cm sized planetesimals can reach densities at least one order of magnitude larger than 
the 1000 cm case. Right panel: corresponding gas density distribution at the end of the simulation, scaled down by two orders of 
magnitude for a direct comparison with the planetesimal density. 



disc surface density does not influence the results of the sim- 
ulations (the planetesimal SPH particles are just a "tracer" 
of the evolution of the solid bodies in the gas disc). How- 
ever, in order to present illustration values in the analysis 
of our results, we will assume that the initial ratio of the 
planetesimal to gas surface densities is 0.01 in all cases. 

Initially all the planetesimals are located in the z = 
plane. However, during the simulation, the random motions 
induced by the gravitational instabilities rapidly stir the 
planetesimal disc up, so that eventually it acquires a finite 
thickness Hp slightly smaller than the gas disc thickness H. 



4 RESULTS 

Fig. El shows the surface density structure of the planetesi- 
mal discs, one outer orbital time after the introduction of the 
planetesimals in the simulations (i.e., after ~ 125 yrs). At 
this stage, most of the disc has evolved for several dynamical 
timescales, so that any initial transient features have disap- 
peared. The figure refers to the cases where the planetesimal 
sizes were 50 cm (left panels) and 1000 cm (right panels). 
The upper panels show the logarithm of the surface density 
Ep (in order to have a direct comparison with Fig. |21 the 
surface densities have been multiplied by 100 and the same 
colour scale has been used). The bottom panels show the 
ratio Sp/E of the planetesimal and gas surface densities. 
The colour scales in the latter plots are exactly the same 
for the two different planetesimal sizes and covers the range 
0.004 < Ep/E < 0.04. 

These plots clearly show how the planetesimal evolu- 
tion changes with changing planetesimal size. As expected, 
the 50 cm planetesimals are strongly influenced by the gas 



drag and display a spiral pattern with very thin spiral arms, 
indicating that the planetesimals are concentrated at the 
bottom of the potential. The effect is reduced in the 1000 
cm case, where the spiral structure in the planetesimal disc 
is similar to that of the gas disc, indicating that the plan- 
etesimals are pushed into the spirals mainly because of the 
gravitational field (note that since the planetesimals have no 
pressure support, we expect the spiral arms to be slightly 
thinner even if no drag force is introduced). A similar struc- 
ture, with relatively broad spiral arms, was indeed also seen 
in the simulation with no drag force. 

The bottom panels in Fig. El show how the concentra- 
tion of planetesimals is modified by the combined effect of 
gas drag and gravity. If the gas and planetesimals respond 
in the same way to the gravitational instabilities, the den- 
sity ratio would be uniform through the disc, so that any 
non-uniformity in the bottom panels is a direct measure of 
the concentration effect on the planetesimals caused by the 
combination of gas gravity and drag force. Clearly, the 50 cm 
planetesimal reach a much higher concentration compared 
to the 1000 cm case. At the end of the simulation Ep/E 
has increased by a factor ~ 3 for the 1000 cm case, while 
in the 50 cm case the maximum increase can be as high as 
~ 50, in some regions therefore reaching surface densities 
comparable to that of the gas. 

Fig-Elshows a histogram of the distribution of planetes- 
imal volume densities, Pp, for both the 50 cm and 1000 cm 
particles (left panel), and a histogram of the distribution of 
gas volume densities, p, scaled down by two orders of mag- 
nitude for a direct comparison with the planetesimal density 
(right panel), at the end of the simulations. In both panels 
A'tot is the total number of particles of the type being consid- 
ered. The planetesimal volume densities were determined by 
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Figure 5. Distribution of volume density ratios (planetesi- 
mal/gas) for the two different sizes considered (solid line: 50 cm; 
dotted line: 1000 cm) at the end of the simulation. While the 
planetesimal/gas density ratio is generally smaller than 0.05 for 
the 1000 cm case (i.e., the concentration enhancement is smaller 
than 5), in the 50 cm case it can reach values of the order of unity 
(concentration enhancement ~ 100). 



assuming, as mentioned earlier, that the initial planetesimal 
to gas surface density ratio was 0.01. The distribution for the 
50 cm planetesimals (solid line - left panel) shows a tail at 
high densities extending to more than an order of magnitude 
above the 1000 cm planetesimals (dotted line - left panel). 
The distribution of gas volume densities (right panel) is also 
similar to that of the 1000 cm particles, again illustrating 
that the 1000 cm particles are more strongly influenced by 
the gravitational potential than by the drag force. The con- 
centration of planetesimals is further illustrated in Fig. |3 
which shows the distribution of the ratio of the planetesi- 
mal volume density to the gas volume density. This shows 
that the volume density of the 50 cm planetesimals (solid 
line) may increase to a value similar to the gas density. The 
volume density ratio for the 1000 cm particles (dotted line), 
on the other hand, barely exceeds 0.1 and in most regions 
has a value below 0.05. Given an "unperturbed" density ra- 
tio of 0.01, this means that the volume density of the 1000 
cm particles is almost never enhanced by more than a fac- 
tor of 10. These results confirm our expectations that for 
the planetesimals with the largest expected radial drift, the 
combined effect of gas drag and spiral structure induced by 
self-gravity leads to a significant concentration of the plan- 
etesimals along the spiral arms. 

Fig. |S| shows how the volume density, as seen by a rep- 
resentative planetesimal located at ~ 8 au, varies with time. 
The left and the right panels display the results for the 50 cm 
and the 1000 cm cases, respectively. The solid line shows the 
planetesimal volume density, while the dotted line shows the 
corresponding gas volume density. The differences between 
the two cases are striking also in this case. The planetes- 
imal density for the 1000 cm particles essentially follows 



closely the gas density, with relatively small variations, and 
oscillates between high and low values as the planetesimal 
goes in and out of the regions of enhanced gas density (the 
spiral arms). In the 50 cm case, the planetesimal density 
reaches extremely high values becoming comparable to the 
gas density when the planetesimal moves into a spiral arm. 
The planetesimal density can also remain high for as long 
as ~ 20 yrs, comparable to the dynamical timescale at 8 au. 

To summarize, Fig.|S]shows that at a given time a sig- 
nificant fraction of the 50 cm planetesimals have a large 
concentration (say, Pp/p > 0.1). On the other hand. Fig. 
|S] shows that a given particle spends only a fraction of the 
total simulation at high Pp/p. This suggest that the total 
fraction of planetesimals which at some stage during the run 
have a large Pp/p is actually larger than the corresponding 
fraction taken at a given time. To estimate this effect we 
have computed the maximum value of Pp/p attained by the 
planetesimals during the whole simulation (in order to re- 
duce computational time, we have performed this analysis 
only for a subset of the total number of planetesimal SPH 
particles). The dotted line in Fig. |7| shows the cumulative 
distribution of Pp/p at a given time (i.e., the cumulative 
distribution corresponding to the solid line of Fig. |^, while 
the solid line shows the cumulative distribution of the maxi- 
mum Pp/p, computed as described above. This figure clearly 
shows that, while at any given time during the simulation 
the fraction of planetesimals with Pp/p > 0.1 (i.e., Pp/p en- 
hanced by more than a factor 10) is ~ 25%, the fraction of 
planetesimal that at some stage during the simulation at- 
tains the same value of Pp/p is significantly higher, ~ 75%. 



5 DISCUSSION AND CONCLUSIONS 

The concentration of planetesimal, due to the combined ef- 
fect of gas drag and gravity, that we find in our simulations 
can have a significant effect on the process of coagulation 
of planetesimal into larger bodies. This can occur either by 
increasing the planetesimal collision rate and/or by making 
the planetesimal sub-disc gravitationally unstable. 

The collision rate of planetesimals is proportional to 
Tipfjp, where rip = pp/rrip is the planetesimal number den- 
sity and (jp is their velocity dispersion. Since in our case the 
Safronov number Q = Gmp/2apa is always much smaller 
than unity, we can neglect the effect of gravitational focus- 
ing. In order to assess the effect of gas drag on the collision 
rate, we have first computed the azimuthally averaged value 
of TipCTp from the simulation with no drag force, (npCTp)nd, as 
a function of radius. We have then computed, for every plan- 
etesimal SPH particle in both the 50 cm case and the 1000 
cm case, the ratio np(jp/(npap)nd, where the average value 
is computed at the same radial location in the disc. If the 
gas drag had no effect on the collision rate, the distribution 
of np(Tp/(npO"p)nd would be strongly peaked around unity. 
Fig-IHlshows the distribution of npcrp/(npcrp)nd that we have 
obtained in the three cases (no drag force: dashed line; 1000 
cm size: dotted line; 50 cm size: solid line). As expected, the 
distribution for the no-drag simulation is strongly peaked 
around unity. Fluctuations of UpGp not related to the gas 
drag result only in an increase of the collision rate by no 
more than a factor ^ 6. In contrast, the introduction of the 
gas drag leads to a broader distribution of collision rates. 
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Figure 6. Time evolution of the densities for a 50 cm planetesimal (left panel) and for a 1000 cm planetesimal (right panel). The solid 
line shows the planetesimal density, while the dashed line shows the gas density. 
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Figure 7. Solid line: cumulative distribution of the maxiiiiuiii 
of the density ratio pp/p attained during the run. Dotted line: 
cumulative distribution of pp/p at a given time. Both plots refer 
to the 50 cm case. While at any given time during the simulation 
the fraction of planetesimals with pp/p > 0.1 (i.e., Pp/p enhanced 
by more than a factor 10) is ~ 25%, the fraction of planetesimal 
that at some stage during the simulation attains the same value 
of Pp/p is significantly higher, ~ 75%. 



Figure 8. Distribution of ripcrp / (ripcrp) for the 50 cm case 
(solid line), for the 1000 cm case (dashed line) and for the no 
drag simulation (dotted line). Fluctuations of ripCTp not related 
to the gas drag result only in an increase of the collision rate by 
no more than a factor ~ 6. In contrast, the introduction of the gas 
drag leads to a broader distribution of collision rates, especially 
for the 50 cm case, where the distribution has a tail extending up 
to two orders of magnitude above the average value. 



especially for the 50 cm case where the distribution has a 
tail extending more than two orders of magnitude above the 
average value. As discussed earlier, the fraction of particles 
concentrated in the spiral arms at a given time is smaller 
than the fraction of particles that, during the course of the 
whole simulation, are at some stage concentrated in the spi- 



ral arms. Since the enhancement in collision rate is due to 
the enhanced density resulting from the concentration of the 
planetesimals in the spiral arms, the number of particles over 
the entire simulation time that at some stage are in a region 
of enhanced collision rate will also be greater than the num- 
ber at a single time. Depending on how well particles of this 
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size stick together during collisions, this enhanced collision 
rate could play an important role in the growth of larger 
planetesimals. 

A planetesimal surface density enhancement of a fac- 
tor ~ 20 may also be sufficient to make the planetesima l 
sub-disc gravitationally unstable (e.g-. lYoudin fc Shul2002f) . 
This is exactly the range of concentrations that we achieve 
in 50 cm simulation. However, since we have neglected the 
planetesimal self-gravity, we are not able to obtain a gravita- 
tional instability in the planetesimal disc in our simulations 
(for a more detailed study of t he stability of a two com - 
ponent self-gravitating disc, see iBertin &: Romeol il988l) '). 
In order to assess the importance of the planetesimal self- 
gravity we have performed two separate tests. 

We have first computed the quantity 

Pp _ ppR'' _ Gpp , , 

where R is the radius of a region that has a local planetesi- 
mal density of pp . This quantity is a measure of the relative 
effects of local gravitational collapse for the planetesimals 
versus tidal disruption. Fig.|^shows the distribution of pp/p 
in the 50 cm case (solid line) and in the 1000 cm case (dot- 
ted line). In the 50 cm case Pp/p reaches much higher values 
than in the 1000 cm case, becoming as large as ~ 1. This sug- 
gest that the density enhancements in the 50 cm simulation 
may be gravitationally significant and that the planetesimal 
disc could become gravitationally unstable. 

As a separate measure, we have also computed the Jeans 
mass, 

at the location of every planetesimal particle. If, at any lo- 
cation, the minimum resolvable mass (the mass, M, con- 
tained within 2 smoothing lengths) exceeds the local Jeans 
mass, then the planetesimals' self-gravity has to play a role. 
The distribution of M/Mj is shown in the right panel of 
Fig- El In the 50 cm case, the distribution is shifted towards 
higher values and is much broader than in the 1000 cm case. 
For the most massive clumps, the ratio M/Mj can become 
comparable to or even larger than unity, indicating that the 
gravitational instability would play a significant role in these 
regions. 

Note that since the vertical scale height of the planetes- 
imal disc is comparable to the scaleheight of the gas disc, the 
high densities achieved in our simulation using 50 cm parti- 
cles is only due to radial and azimuthal compressions rather 
than by the vertical settling of the planetesimals in the mid- 
plane. Therefore, u nlike in the standard picture for gravi- 
tatio nal instability jGoldreich fc WarJl973t lYoudin fc Shi 
l2002tl . our results should not be affected by additional tur- 
bulence generated by the vertical shear between the gas and 
the dust. 

To summarize, in this work we have shown how the in- 
teraction between a self-gravitating gaseous protoplanetary 
disc and embedded planetesimals plays an important role in 
accelerating planetesimal growth. However, there are a num- 
ber of important effects that we have neglected in this first 
approach to the problem, (i) The planetesimals have been 
essentially modelled as test particles, ignoring both the back 
reaction of the planetesimals on the gas and the planetesi- 



mal self-gravity. In regions where the planetesimal density is 
enhanced (becoming in some regions comparable to the gas 
density) both effects may be very important, (ii) In each 
simulation the planetesimals are assumed to be of a single 
size and the volume density is computed by assuming that 
the initial surface density ratio is 0.01. In reality there will 
be a range of plan etesimal sizes (e.g., iMathis et al.lll97^ : 
iMizuno et alj [l98^: note, however, that these studies con- 
sider grain sizes much smaller than those considered here) 
and only those with sizes, in this case, between ~ 10 and 100 
cm (see Fig. will be significantly infiuenced by the self- 
gravitating structures in the gas disc. During the evolution 
of the protoplanetary nebula there might well be some stage 
where most of the planetesimal mass is contained within 
a relatively small range of sizes, so that our assumptions 
may not be unrealistic. When this size range includes the 
size for which the drift induced by the drag is significant, 
we can expect a significant increase in collision rate and 
an enhanced tendency toward gravitational collapse. To ad- 
dress the details of these processes we would need to consider 
many other effects (such as the sticking properties of plan- 
etesimals) which are beyond the scope of the present paper. 
(in) In this work we have considered, as an illustration, the 
self-gravitating structure resulting from a relatively massive 
gas disc (with a mass ~ 0.25A/j,). However, a well defined 
spiral structure is also present in significantly less massive 
discs (Lodato fc Rice 200^. It would then be interesting to 
check the dependence of some of the details of the results 
described here (such as what is the relevant size range for 
the planetesimals that display the strongest response to the 
spiral structure) on the specific choice of the gas disc prop- 
erties. 

We plan in the future to include some of the effects 
described above, such as the backreaction of the planetesi- 
mals on the disc gas, the effect of planetesimal self-gravity, 
and we plan to consider various disc masses. However, it 
seems clear that if protoplanetary discs experience a self- 
gravitating phase, the resulting disc structures could well 
play an important role in planetesimal evolution and growth 
and could ultimately influence the growth of terrestrial plan- 
ets and the cores of gas/ice giant planets. Also, since a pro- 
toplanetary disc is most likely to become gravitationally un- 
stable early in the star formation process, we might expect 
substantial processing of the dust prior to the optically vis- 
ible Classical T Tauri phase. 
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Figure 9. Left panel: distribution of pp/p for the 50 cm (solid line) and the 1000 cm case (dotted line). Right panel: distribution of 
M/Mj in the two cases. These results show that some of the high density clumps observed in the 50 cm case should be subject to 
gravitational instability. 
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